#include"../header/gaussLegendre.h"
#define _output_

void gaussLegendreSeq(FILE * arquivo) {

#ifdef _output_
	int i = 0;

	printf("Inicializando Gauss-Legendre...\n");
#endif

	mpf_t A0, B0, T0, P0, AN, BN, TN, PN, Pi;
	mpf_t diferenca, precisao;
	mpf_t aux;

	struct timeval inicial, final;

	int bits = 33219281; /* 10^7 * log(10) */ 

	mpf_set_default_prec(bits);

	mpf_init(A0); 
	mpf_init(B0);
	mpf_init(T0);
	mpf_init(P0);
	mpf_init(AN);
	mpf_init(BN);
	mpf_init(TN);
	mpf_init(PN);

	mpf_init(Pi);
	mpf_init(diferenca);
	mpf_init(precisao);

	mpf_init(aux);

	inicial = tempo();
	
	/* A0 = 1 */
	mpf_set_d(A0, 1.0f);

	/* B0 = 1/sqrt(2) */
	mpf_set_d(B0, 1.0f);
	mpf_set_d(aux, 2.0f);
	mpf_sqrt(aux, aux);
	mpf_div(B0, B0, aux);

	/* T0 = 1/4 */
	mpf_set_d(T0, 0.25f);

	/* P0 = 1 */
	mpf_set_d(P0, 1.0f);

	/* Precisao = 1/10^(10^7) */
	mpf_set_d(precisao, 10);
	mpf_pow_ui(precisao, precisao, (unsigned long int)10000000);
	mpf_ui_div(precisao, (unsigned long int) 1, precisao);

#ifdef _output_
	printf("Inicializado.\n");
#endif

	do {

#ifdef _output_
		printf("Iteracao %d.\n", i);
		i++;
#endif

		/* AN = (A0 + B0)/2 */
		mpf_add(aux, A0, B0);
		mpf_div_ui(AN, aux, (unsigned long int) 2);

		/* BN = sqrt(A0 * B0) */
		mpf_mul(aux, A0, B0);
		mpf_sqrt(BN, aux);		

		/* TN = T0 - P0 * (A0 - AN)^2 */
		mpf_sub(aux, A0, AN);
		mpf_pow_ui(aux, aux, (unsigned long int) 2);
		mpf_mul(aux, aux, P0);
		mpf_sub(TN, T0, aux);

		/* PN = 2 * P0 */
		mpf_mul_ui(PN, P0, (unsigned long int) 2);

		/* Prepara proxima iteracao.  */
		mpf_set(A0, AN);
		mpf_set(B0, BN);
		mpf_set(T0, TN);
		mpf_set(P0, PN);

		/* Calcula diferenca. (diferenca = AN - BN) */
		mpf_sub(diferenca, AN, BN);
		mpf_abs(diferenca,diferenca);


	} while(mpf_cmp(diferenca, precisao) > 0); /* Precisao atingida! */	

	/* PI = (AN + BN)^2 / (4 * TN) */
	mpf_mul_ui(aux, TN, 4.0f);
	mpf_add(Pi, AN, BN);
	mpf_pow_ui(Pi, Pi, (unsigned long int) 2);
	mpf_div(Pi, Pi, aux);

	final = tempo();

    printf("Tempo de execucao: %.3lf s\n", difftempo(final,inicial));

	/* Imprime PI... */
	mpf_out_str(arquivo, 10, 10000001, Pi); 

	
	mpf_clear(AN);
	mpf_clear(A0);
	mpf_clear(BN);
	mpf_clear(B0);
	mpf_clear(PN);
	mpf_clear(P0);
	mpf_clear(TN);
	mpf_clear(T0);
	mpf_clear(Pi);

	mpf_clear(diferenca);
	mpf_clear(precisao);
	mpf_clear(aux);
	
}

mpf_t A0_par, A1_par, B0_par, T0_par, P0_par, AN_par, BN_par, TN_par, PN_par, Pi_par;
mpf_t diferenca_par, precisao_par;
int convergiu_gauss;

void * calc_a_gauss()
{
	/* AN = (A0 + B0)/2 */
	mpf_add(AN_par, A1_par, B0_par);
	mpf_div_ui(AN_par, AN_par, (unsigned long int) 2);

	return NULL;
}

void * calc_b()
{
	/* BN = sqrt(A0 * B0) */
	mpf_mul(BN_par, A1_par, B0_par);
	mpf_sqrt(BN_par, BN_par);	
	
	return NULL;
}

void * calc_t()
{
	/* TN = T0 - P0 * (A0 - AN)^2 */
	mpf_sub(TN_par, A0_par, A1_par);
	mpf_pow_ui(TN_par, TN_par, (unsigned long int) 2);
	mpf_mul(TN_par, TN_par, P0_par);
	mpf_sub(TN_par, T0_par, TN_par);

	return NULL;
}

void * calc_p()
{
	/* PN = 2 * P0 */
	mpf_mul_ui(PN_par, P0_par, (unsigned long int) 2);

	return NULL;
}

void * checa_pi()
{
	/* Calcula diferenca. (diferenca = AN - BN) */
	mpf_sub(diferenca_par, A1_par, B0_par);
	mpf_abs(diferenca_par,diferenca_par);

	if( mpf_cmp(diferenca_par, precisao_par) < 0 )
		convergiu_gauss = 1;

	return NULL;
}

void gaussLegendrePar(FILE * arquivo)
{
	printf("Inicializando Gauss-Legendre...\n");

	int i = 0;

	struct timeval inicial, final;

	pthread_t a_thread, b_thread, t_thread, p_thread, pi_thread;	

	int bits = 33219281; /* 10^7 * log(10) */ 

	mpf_t aux;

	mpf_set_default_prec(bits);

	convergiu_gauss = 0;

	mpf_init(A0_par);
	mpf_init(A1_par); 
	mpf_init(B0_par);
	mpf_init(T0_par);
	mpf_init(P0_par);
	mpf_init(AN_par);
	mpf_init(BN_par);
	mpf_init(TN_par);
	mpf_init(PN_par);

	mpf_init(Pi_par);
	mpf_init(diferenca_par);
	mpf_init(precisao_par);

	mpf_init(aux);

	inicial = tempo();

	/* A0 = 1 */
	mpf_set_d(A0_par, 1.0f);

	/* B0 = 1/sqrt(2) */
	mpf_set_d(B0_par, 1.0f);
	mpf_set_d(aux, 2.0f);
	mpf_sqrt(aux, aux);
	mpf_div(B0_par, B0_par, aux);

	/* T0 = 1/4 */
	mpf_set_d(T0_par, 0.25f);

	/* P0 = 1 */
	mpf_set_d(P0_par, 1.0f);

	/* Precisao = 1/10^(10^7) */
	mpf_set_d(precisao_par, 10);
	mpf_pow_ui(precisao_par, precisao_par, (unsigned long int)10000000);
	mpf_ui_div(precisao_par, (unsigned long int) 1, precisao_par);
	
	printf("Inicializado.\n");

	/*primeira iteracao*/
	mpf_set(A1_par,A0_par);

	pthread_create(&a_thread,NULL,calc_a_gauss,NULL);
	pthread_create(&b_thread,NULL,calc_b,NULL);

	pthread_join(a_thread,NULL);
	pthread_join(b_thread,NULL);

	mpf_set(B0_par,BN_par);
	mpf_set(A1_par, AN_par);

	/*segunda iteracao*/
	pthread_create(&a_thread,NULL,calc_a_gauss,NULL);
	pthread_create(&b_thread,NULL,calc_b,NULL);
	pthread_create(&t_thread,NULL,calc_t,NULL);
	pthread_create(&p_thread,NULL,calc_p,NULL);

	pthread_join(a_thread,NULL);
	pthread_join(b_thread,NULL);
	pthread_join(t_thread,NULL);
	pthread_join(p_thread,NULL);

	mpf_set(A0_par, A1_par);
	mpf_set(A1_par, AN_par);
	mpf_set(B0_par, BN_par);
	mpf_set(T0_par, TN_par);
	mpf_set(P0_par, PN_par);

	
	/*pipeline*/	

	do {

		printf("Iteracao %d.\n", i);
		i++;

		pthread_create(&a_thread,NULL,calc_a_gauss,NULL);
		pthread_create(&b_thread,NULL,calc_b,NULL);
		pthread_create(&t_thread,NULL,calc_t,NULL);
		pthread_create(&p_thread,NULL,calc_p,NULL);
		pthread_create(&pi_thread,NULL,checa_pi,NULL);

		pthread_join(a_thread,NULL);
		pthread_join(b_thread,NULL);
		pthread_join(t_thread,NULL);
		pthread_join(p_thread,NULL);
		pthread_join(pi_thread,NULL);

		/* Prepara proxima iteracao.  */
		mpf_set(A0_par, A1_par);
		mpf_set(A1_par, AN_par);
		mpf_set(B0_par, BN_par);
		mpf_set(T0_par, TN_par);
		mpf_set(P0_par, PN_par);

	} while(!convergiu_gauss); /* Precisao atingida! */	

	/* PI = (AN + BN)^2 / (4 * TN) */
	mpf_mul_ui(aux, TN_par, 4.0f);
	mpf_add(Pi_par, A1_par, B0_par);
	mpf_pow_ui(Pi_par, Pi_par, (unsigned long int) 2);
	mpf_div(Pi_par, Pi_par, aux);
	
	final = tempo();

    printf("Tempo de execucao: %.3lf s\n", difftempo(final,inicial));

	/* Imprime PI... */
	mpf_out_str(arquivo, 10, 10000001, Pi_par); 
	
	mpf_clear(AN_par);
	mpf_clear(A0_par);
	mpf_clear(A1_par);
	mpf_clear(BN_par);
	mpf_clear(B0_par);
	mpf_clear(PN_par);
	mpf_clear(P0_par);
	mpf_clear(TN_par);
	mpf_clear(T0_par);
	mpf_clear(Pi_par);

	mpf_clear(diferenca_par);
	mpf_clear(precisao_par);
	mpf_clear(aux);

}

